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Abstract. Recent work has explored solver strategies for the linear system of equations arising 
from a spectral Galerkin approximation of the solution of PDEs with parameterized (or stochastic) 
inputs. We consider the related problem of a matrix equation whose matrix and right hand side 
depend on a set of parameters (e.g. a PDE with stochastic inputs semidiscretized in space) and 
^s^j ■ examine the linear system arising from a similar Galerkin approximation of the solution. We derive 

a useful factorization of this system of equations, which yields bounds on the eigenvalues, clues to 
preconditioning, and a flexible implementation method for a wide array of problems. We complement 
this analysis with (i) a numerical study of preconditioners on a standard elliptic PDE test problem 
and (ii) a fluids application using existing CFD codes; the MATLAB codes used in the numerical 
studies arc available online. 
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1. Introduction. Complex engineering models are often described by systems of 
equations, where the model outputs depend on a set of input parameters. Given values 
for the input parameters, the model output may be computed by solving a discretized 
differential equation, which typically involves the solution (or series of solutions) of a 
system of linear equations for the unknowns. Each of these computations may be very 
expensive depending on factors such as grid resolution or physics components in the 
model. As the role of simulation gains prominence in areas such as decision support, 
design optimization, and predictive science, understanding the effects of variability 
^ 1 in the input parameters on variability in the model output becomes more important. 

Exhaustive parameter studies (e.g. uncertainty quantification, sensitivity analysis, 
model calibration) can be prohibitively expensive - particularly for a large number of 
input parameters - and therefore accurate interpolation and robust surrogate models 
are essential. 

We consider the model problem of a matrix equation whose matrix and right hand 
side depend on a set of parameters. Let s e [— 1, l] d be a set of input parameters 
from an idealized parameter space. We equip the parameter space with a bounded, 
separable, positive weight function uj : [—1, l] d i->- R + , where tu(s) — u>i(si) ■ ■ -ujd(sd)- 
(In a probabilistic context, this weight function represents a probability measure on 
the input parameter space.) Let A(s) be an N x N matrix-valued function where we 
assume that A(s) is invertible for all s E [— and let b(s) be the vector-valued 

function with N square-integrable components. We seek the vector valued function 
x(s) that satisfies 

A(s)x(s) = b(s), se[-l,l] d . (1.1) 

Such parameterized matrix problems often arise as an intermediate step when com- 
puting an approximate solution of a complex model with multiple input parameters. 
They appear in such diverse fields as differential equations with random (or parame- 
terized) inputs [U [31] , electronic circuit design [23] , image deblurring models [5] , and 
ranking methods for nodes in a graph [25[[7]. Given a parameterized matrix equation, 
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one may wish to approximate the vector valued function that satisfies the equation 
or estimate its statistics. Once an approximation is constructed that is cheaper to 
evaluate than the true solution, statistics - such as mean and variance - of the ap- 
proximation represent estimates of statistics of the true solution. 

The approximation model is a vector of multivariate polynomials represented as 
a series of orthonormal polynomial basis functions where each basis function is a 
product of univariate orthonormal polynomials. We employ the standard multi-index 
notation; let a — (ai, . . . , ay) G N d be a multi-index, and define the basis polynomial 

7r a (s) =7r ai (s 1 )---7r a< j(s d ). (1.2) 

The polynomial ir ai (si) is the orthonormal polynomial of degree ai, where the or- 
thogonality is defined with respect to the weight function uji(si). Then for a, j3 £ N d , 

J hl l]d ^WMsMs) ds = { ^ " t =f w . se (1.3) 

where equality between multi-indicies means component-wise equality. For a given 
index set X C N d with size \L\ < oo, the polynomial approximation can be written 

x ( s ) ~ ^2x a Tr a (s) = Xtt(s). (1.4) 

The N- vector x Q is the coefficient of the series corresponding to 7r Q (s). The N x \1\ 
matrix X has columns x a , and the parameterized vector 7r(s) contains the basis 
polynomials. The goal of the approximation method is to compute the unknown 
coefficients X. 

Such polynomial models have become popular for approximating the solution 
of PDEs with random inputs; they appear under names such as polynomial chaos 
methods [3S] , stochastic finite element methods [TB] , stochastic Galerkin methods [2] , 
and stochastic collocation methods [35] [3] . The Galerkin methods compute the series 
coefficients such that the equation residual is orthogonal to the approximation space 
defined by the index set I; they typically employ a full polynomial basis of order n 
given by 

X = l n = {a e N d : oH h a d < n}, (1.5) 

although such a basis set is not strictly necessary. The number of terms in this basis 
set is \I n | = ( n ^ d ) , which grows exponentially with the number of input parameters d. 
The process of computing the coefficients involves solving a linear system of size N\I\ x 
N\X\ , which can be prohibitively large for even a moderate number of input parameters 
(6 to 10) and low order polynomials (degree > 5). For this reason, there has been 
a flurry of recent work on solver strategies for the matrix equations arising from the 
Galerkin methods [27l [32l [211 EH H21 Ej , including papers on preconditioning |30l El 
[29l [33] . Such work has relied on knowing the matrix valued coefficients A a of a series 
expansion of the parameterized matrix A(s), 

A(s) =^A Q ^ Q ( S ), (1.6) 

a 

which typically come from a specific form of the coefficients of the elliptic operator in 
the PDE models. Another drawback of the Galerkin method is its limited ability to 
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take advantage of existing solvers for the problem A(X)x(X) = 6(A) given a parame- 
ter point A e [—1, l] d . In contrast, the pseudospectral and collocation methods can 
compute coefficients of the model (|1.4p using only evaluations of the solution vector 
x(X). The chosen points typically correspond to a multivariate quadrature rule, and 
computing the coefficients of the polynomial model (| 1 .4[) becomes equivalent to ap- 
proximating its Fourier coefficients with the quadrature rule [3j [8] . Thus, from the 
point of view of code reuse and rapid implementation, the pseudospectral/collocation 
methods have a distinct advantage. 

We propose a variant of the Galerkin method that alleviates the drawbacks of 
limited code reuse and memory limitations. By formally replacing the integration in 
the Galerkin method by a multivariate quadrature rule - a step that is more often 
than not performed in practice - we derive a decomposition of the linear system of 
equations used to compute the Galerkin coefficients of (|1.4p ; such a method has been 
called Galerkin with numerical integration in the context of numerical methods for 
PDEs [4 . This decomposition allows us to compute the Galerkin coefficients using 
only evaluations of the parameterized matrix A(X) and parameterized vector 6(A) 
for points A G [—1, l] d corresponding to a quadrature rule. In fact, if one requires 
only matrix-vector multiplies as in Krylov-based iterative methods (e.g. CG [5D] or 
minres [26] solvers) for the Galerkin system of equations, this restriction can be 
relaxed to matrix-vector products A(\)v for a given TV-vector v; there is no need 
for the coefficients of an expansion of A(s) as in (jl.6|) . Therefore the method takes 
full advantage of sparsity of the parameterized system resulting in reduced memory 
requirements. The decomposition also yields insights for preconditioning the Galerkin 
system that generalize existing work. Additionally, the decomposition immediately 
reveals bounds on the eigenvalues of the Galerkin system for the case of symmetric 



The remainder of the paper is structured as follows. In Section [5J we derive the 
Galerkin method for a given problem and basis set. We then derive the decomposition 
of the matrix used to compute the Galerkin coefficients and examine its consequences 
including bounds on the eigenvalues of the matrix, strategies for simple implementa- 
tion and code reuse, and insights into preconditioning. We then provide a numerical 
study of various preconditioners suggested by the decomposition using the common 
test case of an elliptic PDE with parameterized coefficients in Section [3J the codes 
for the numerical study are available in a Matlab suite of tools accessible online. 
To emphasize the advantages of code reuse, we apply the method to an engineering 
test problem using existing codes in Section [4] Finally, we conclude in Section [5] with 
summarizing remarks. 

1.1. Notation. For the remainder of the paper, we will use the bracket notation 
(•) to denote a discrete approximation to the integral with respect to the weight 
function ui. In other words, for functions / : [—1, l] d — > R, 



where the points \p = (A^, . . . ,\p d ) € [—1, l] d and weights vp € R dchne a mul- 
tivariate quadrature rule for multi-indicies in the set J C N d . We will abuse this 
notation by putting matrix and vector valued functions inside the brackets, as well. 
For example, (A) denotes the mean of A(s) computed with a quadrature rule. 



A(s). 




(1.7) 
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2. A Spectral Galerkin Method. We first review the Galerkin method [5] for 
computing the coefficients of the polynomial model (|1.4I) . Define the residual 

r(y,s) = A(s)y(s)-b(s), (2.1) 

and let x g (s) be the Galerkin approximation. We require that each component of the 
residual r.i(x g , s) be orthogonal to the approximation space defined by the span of ir a 
with a £ X, 

(n{x g )n a ) = 0, i = l,...,N, a£X. (2.2) 

We can combine the equations in (|2.2[) using the matrix notation 

(r(x g )n T ) = ((Ax g -b)n T ) = 0, (2.3) 

or equivalently, upon substituting the model x g (s) = X-7r(s), 

(AX.tztt t ) = (bir T ) . (2.4) 

Using the vec notation [T9| Section 4.5], we can rewrite (|2.4[) as 

(tvtv t ® A) x = (tt ® b) . (2.5) 

where x = vec(X) is an iV|X| x 1 constant vector equal to the columns of X stacked 
on top of each other. The constant matrix (tv-k t <g) A} has size iV|I| x N\X\ and a 
distinct block structure; the a, /3 block of size N x N is equal to (n a TTpA) for multi- 
indicies a, (3 £ X. Similarly, the a block of the iV|I| x 1 vector (ir ® b) is equal to 
(bwa), i.e. the Fourier coefficient of b associated with n a (s) approximated with the 
quadrature rule. 

Much of the literature on PDEs with random inputs [23 [30] points out an in- 
teresting block sparsity pattern that arises in the matrix (mr T ® A) when the pa- 
rameterized matrix A(s) depends at most linearly on any componenets of s; this is 
a result the orthogonality of the bases 7r(s). For general analytic dependence on the 
parameters, such sparsity patterns do not appear [14]. However, one can always mim- 
ick the sparsity pattern of A(s) with a simple reordering of the variables. By taking 
the transpose of (|2.4|) and using the same vec operations, if we define x = vec(X T ) 
(i.e. the same unknowns reordered), then 

(A (g) 7T7r T ) x = (6(g) tt) . (2.6) 

Notice that the matrix in (|2.6I) retains the sparsity of A(s) in its blocks, since each 
i,j block of size \X\ x \X\ is equal to (Aijinr T }, where Aij(s) is the i,j element of 
A(s). For the remainder of the paper, however, we will work with the form (|2.5[) . 

By writing out the numerical integration rule for the integrals used to form 
(tttt t ® A), we uncover an interesting decomposition, which we state as a theorem. 

Theorem 2.1. Let {{Xp,vp)~} with j3 £ J be a multivariate quadrature rule. The 
matrix (irir T ® A} can be decomposed as 

(mr T ® A) = (Q(g)I)A(A)(Q(gI) T , (2.7) 

where I is the N x N identity matrix, and Q is a matrix of size \X\ x \ J\ - one row 
for each polynomial basis and one column for each point in the quadrature rule. The 
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matrix A(X) is a block diagonal matrix of size N\J\ x iV|J7"| where each nonzero block 
is A(Xp) for /3 G J. 

Proof. Writing out the quadrature rules, 

(ttt^ ® A) = [ir(Xfi)ir(X fi ) T ® A(\p)] Pp. (2.8) 

Notice that the elements of the vector ir(Xp) are the polynomials Tr a (s) evaluated at 
the quadrature point Xp. If we define the vectors 

= y/UpiriXp), (2.9) 

then we have 

(7T7r T <g> A) = Y ®A{Xp). (2.10) 

Let Q be the matrix whose columns are q^, and define the block diagonal matrix 
A(X) with diagonal blocks equal to A(Xp). Then for an TV x TV identity matrix I, we 
can rewrite (|2.10|) as (|2.7|) . as required. □ 

As an aside, we note that if A(s) depends polynomially on the parameters s, 
then each integrand in the matrix {t^tt t ® A) is a polynomial in s. Therefore by the 
polynomial exactness, there is a Gaussian quadrature rule such that the numerical 
integration approach exactly recovers the true Galerkin matrix. 

The elements of Q are the orthogonal polynomials evaluated at the quadrature 
points and multiplied by the square root of the quadrature weights. Thus, for two 
polynomials 7r Q (s) and ftp(s) with a,/3el, define r a and rp to be the corresponding 
rows of Q; then 

TaY 1 = ^2 7Ta(A 7 )7T£(A 7 )i/ 7 . (2.11) 

If the quadrature rule is a tensor product Gaussian quadrature rule of sufficiently 
high order to exactly compute the polynomial integrand, then this implies QQ T = I, 
where I is the \X\ x \X\ identity matrix; we will assume from here on that the chosen 
quadrature rule yields this property. We will also assume that the number of basis 
polynomials is less than the number of points used in the quadrature rule, i.e. \X\ < 
\J~\; otherwise the Galerkin matrix (ttit t ® A) will be rank deficient. 

Theorem 12.11 reveals bounds on the eigenvalues of the matrix (tttt t ® A) for the 
case when A(s) is symmetric; we state this as a corollary. 

COROLLARY 2.2. Suppose A(s) is symmetric for all s G [—1, l] d . The eigenvalues 
of (iztt t (g> A) satisfy the bounds 



mm 



(A(Xp))] < 0((7T7r T ®A)) < max[e max (A(A^))] , (2.12) 



where 9(X) denotes the eigenvalues of a matrix X , and # m in(^0 and max {X) denote 
the smallest and largest eigenvalues of X , respectively. 

Proof. Let {(A^, vp)} be a tensor product Gaussian quadrature rule of sufficiently 
high order so that QQ T = I. Let 7r(s) be a vector containing the polynomial bases 
such that [k(s) t , 7r(s) T ] T contains the tensor product polynomial basis corresponding 
to the tensor product quadrature grid {A^}. Let Q be the matrix with columns 
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q_p — Tr(\p)^vp] by construction, the rows of Q are orthonormal and QQ = 0. 
Also, by the Christoffcl-Darboux formula Q21 Theorem 1.32], Q T Q + Q T Q = I. 
Then the matrix 



Z = 



(2.13) 



is square and orthogonal, i.e. Z _1 = Z T . Notice that the matrix (Q® I)A(A)(Q<g)I) T 
is the first \I\ x \I\ principal minor of (Z ® I)A(\)(Z eg) I) T . Then by the interlacing 
theorem [121 Theorem 8.1.7], the eigenvalues of (Q <X> I)A(A)(Q ® I) T are bounded 
by the extreme eigenvalues of (Z ® I)A(\)(Z ® I) T . But since Z T = Z _1 , this is a 
similarity transformation with ^4(A), which completes the proof. □ 

2.1. Iterative Solvers. Due to their size and sparsity, the preferred way of 
solving (|2.5p is with Krylov based iterative solvers [27] that rely on matrix- vector 
products with the matrix (tvk t <£) A}. By employing the decomposition in Theorem 
12.11 we can compute these using only multiplication of the parameterized matrix 
evaluated at the quadrature point by a given vector. More precisely, given a vector 
u = vcc(U), suppose we wish to compute 

v = vcc(V) = (Q<g>I)A(A)(Q<g) I) T u. (2.14) 



We accomplish this in three steps using the properties of the Kronecker product: 

1. W = UQ. Let be a column of W with j3 e J '. 

2. For each j3, yp — A(Xp)wp. Define Y to be the matrix with columns yp. 

3. V = YQ T . 

Step Q] can be thought of as pre-processing, and step [3] as post-processing. In practice, 
each row of Q may have a Kronecker structure corresponding to the tensor product 
quadrature rule. In this case, steps [1] and [3] can be computed accurately and efficiently 
using multiplication methods such as [15] . The second step requires only constant 
matrix-vector products where the matrix is A{s) evaluated at the quadrature points. 
Therefore we can take advantage of a memory-efficient, reusable interface for the 
matrix-vector multiplies that will exploit any sparsity in the matrix. We reiterate 
that this can be accomplished without any knowledge of the specific type of parameter 
dependence in A(s). 

2.2. Preconditioning Strategies. The number of iterations required to achieve 
a given convergence criterion (e.g. a sufficiently small residual) can be greatly reduced 
for Krylov-based iterative methods with a proper preconditioner. In general, precon- 
ditioning a system is highly problem dependent and begs for the artful intuition of 
the scientist. However, the structure revealed by the decomposition from Theorem 
12.11 offers a number of useful clues. 

Suppose we have an N x N matrix P that is easily invertible. We can construct a 
block-diagonal preconditioner I<Ei P _1 , where I is the identity matrix of size \I\ x \I\. 
If we premultiply the preconditioner against the factored form of (t^-k t ® A), we get 

(I <g> P" 1 )(Q ® I)A(A)(Q ® I) T = (Q <g> I)(I <g> P~ 1 )A(A)(Q <g> I) T . (2.15) 

Notice that by the mixed product property and commutativity of the identity matrix, 
the block-diagonal preconditioner slips past Q(g)I to act directly on the parameterized 
matrix evaluated at the quadrature points. The blocks on the diagonal of the inner 
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matrix product are P~ 1 A(X/3) for (3 £ J . In other words, we can choose one constant 
matrix P to affect the parameterized system at all quadrature points. 

A reasonable and popular choice is the mean P = (A); see [3S1[57] for a detailed 
analysis of this preconditioner for stochastic finite element systems. Notice that this 
is also the first N x N block of the Galerkin matrix. However, if A(s) is very large 
or has some complicated parametric dependence, then forming the mean system and 
inverting it (or computing partial factors) for the preconditioner may be prohibitively 
expensive. If the dependence of A(s) on s is close to linear, then P = A((s)) may be 
much easier to evaluate and just as effective. 

One goal of the preconditioner is reduce the condition number of the matrix, and 
one way of achieving this is to reduce the spread of the eigenvalues. If we knew a pri- 
ori which region of the parameter space produced the extrema of the parameterized 
eigenvalues of A(s) (e.g. the boundaries), then we could choose an appropriate pa- 
rameter value to construct an effective preconditioner. Unfortunately, we only get to 
use one such evaluation. Therefore, if the largest possible value of the parameterized 
eigenvalues is very large, we may choose this parameter value. Alternatively, if the 
smallest eigenvalue is close to zero (for positive definite systems), then this may be 
a better option to reduce the condition number of the Galerkin system. In the next 
section, we explore a few choices for the preconditioner P on a standard test problem. 

3. Preconditioning Study. Consider the following parameterized elliptic par- 
tial differential equation with homogeneous Dirichlet boundary conditions; variations 
of this problem can be found in [Tl [2j [16] , amongst others. We seek a solution u(x, s) 
that satisfies 

V • (a(x, s)Vu(x, s)) = 1, ie[0,lf, (3.1) 

where u — on the boundary of the square domain [0, l] 2 , and the parameter space 
is the hypercube s G [—1, 1] 4 equipped with a uniform measure. The logarithm of 
the elliptic coefficient is given by a truncated Karhunen-Loeve like expansion |24| of 
a zero-mean random field with covariance 

C(x u x 2 ) = 2exp(-|| 2 : 1 -x 2 || 2 /2), (3.2) 

so that 

4 

log(a(x, s)) = 2 ^ °kipk{x)s k . (3.3) 
fc=i 

We discretize (|3 . 1 [) using the finite element method implemented in Matlab's PDE 
toolbox on an irregular mesh of N = 1,921 triangles. We solve the eigenvalue prob- 
lem with the discrete covariance matrix to compute the values of the eigenfunctions 
ipk{x) on the given mesh. The square roots of first ten eigenvalues o~k are plotted in 
Figure 13.11 we truncate the expansion at d = 4 to retain roughly 90% of the energy 
of the field. Given a point in the parameter space, the PDE Toolbox allows us to 
access the stiffness matrix. Therefore we can apply Matlab's MINRES solver using 
only matrix-vector multiplies against the stiffness matrix evaluated at the quadrature 
points. For the polynomial basis, we use the normalized multivariate product Legen- 
dre polynomials of order 5, which includes \X\ = 126 basis functions. Therefore the 
number of unknowns is N\I\ — 242,046. We use a tensor product Gauss-Legendre 
quadrature rule of order 12 in each parameter (20,736 points) to implicitly form the 
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Fig. 3.1. Eigenvalues of the KL expansion. The numbers indicate the percentage of field energy 
captured when using that many terms in the expansion. 
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Galerkin matrix; this is more than sufficient to maintain the orthogonality in the basis 
functions. 

To test the preconditioning ideas suggested by the decomposition from Theorem 
I2.1l and equation (|2.f 5[) . we try five different choices of P. For each P, we precomputc 
the Cholesky factors to apply (|2 . 15[) efficiently Figure I3.2I and Table 13 . 1 1 summarize 
the results, which we now explain in detail. In the first type of preconditioning, we 
set P = A(s r ) for a random point s r in [— 1, l] 4 . To sample the likely effectiveness 
of any random point, we select 25 random points and evaluate each choice of s r . 
Second, we set P = A(s max ) where -A(s max ) is the matrix with the largest eigenvalue 
in [—1, l] 4 . To compute s max , we use up to 100 iterations of the power method to 
estimate the largest eigenvalue at each point in the quadrature rule, that is for each 
A/3 for /3 € J\ we also evaluate A{s) for the tensor product of the endpoints as well. 
Third, we set P = A(s m - m ) where A(s m ; n ) is the matrix with the smallest eigenvalue in 
[—1, l] 4 . For this computation, we use Matlab's optimization routine fmincon and 
use eigs/ARPACK to evaluate the smallest eigenvalue [22]. Fourth, we set P = A(s m id) 
where s m i<j is the midpoint of the domain, which is the origin for our experiments. 
Fifth, and finally, we set P = (A), the mean preconditioner. To evaluate the mean, we 
use a 2nd order quadrature rule for a fast approximation and a 5th order quadrature 
rule for a more exact approximation. 

In Figure 13.21 we show the convergence of the 2-norm of the residual for each of 
the preconditioning strategies, as well as no preconditioning. We only show the result 
from the 5th order mean based preconditioner as there was no appreciable difference 
in convergence. This plot clearly shows that the mean and midpoint preconditioners 
are excellent choices. Further note that any random point is more than two orders 
of magnitude better than no preconditioning at all. In fact, using a random point is 
better than using the point with largest eigenvalue and may compare with using the 
point with the smallest eigenvalue. 

Next, Table [3~T1 shows a few timing results from these experiments. We present the 
time taken to compute/setup the preconditioner, the average time between iterations, 
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Fig. 3.2. Convergence of the 2-norm of the residual for a variety of preconditioners. 



Table 3.1 

Timing results for preconditioning showing only two of the random points. 



Method 


Setup (s) 


Iterations 


Avg. Iter, (s) 


Total (s) 


no preconditioning 


0.0 


7075 


6.5 


46251.1 


largest eigenvalue 


176.0 


413 


6.6 


2810.7 


smallest eigenvalue 


16.0 


216 


6.5 


1469.4 


random 


0.1 


238 


6.7 


1661.8 


random 


0.9 


216 


8.4 


2030.2 


mean (5th order) 


4.7 


111 


6.7 


807.2 


mean (2nd order) 


0.7 


110 


6.8 


820.4 


midpoint 


0.1 


110 


7.0 


843.2 



and the total time taken by the MINRES method (excluding preconditioncr setup) . The 
times are from Matlab 2010a on an Intel Core i7-960 processor (3.2 GHz, 4 cores) 
with 24GB of RAM. Matlab used four cores for its own multithreading and we used 
the Parallel Computing toolbox's parf or construction to parallelize the matrix- vector 
product over 4 cores as well. By watching a processor performance meter, we observed 
high utilization of all four cores. We note, however, that during some of the random 
point evaluations, another process intruded and caused the system to swap memory 
to disk heavily. We believe this caused the inconsistent average iteration time. From 
these timing results, we conclude that preconditioning changes the iteration time only 
slightly - if at all. Furthermore, the setup times are often small when compared with 
the savings in runtime. 

The Matlab codes for this study include a utility for generating realizations of 
random fields [2] and a suite of tools called PMPack (Parameterized Matrix Pack- 
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Fig. 4.1. Mesh used to discretize the domain O and compute temperature distribution. 

age) for using spectral methods to approximate the solution of parameterized matrix 
equations. The codes for generating the results of the study can be found at: 
http : //www. Stanford. edu/~dgleich/publicat ions/20 10/ codes/pmpack-sisc/ 

4. Application Heat Transfer with Uncertain Material Properties. 

As a proof of concept, we examine an application from computational fluid dynamics 
with uncertain model inputs. The flow solver used to compute the deterministic 
version of this problem - i.e. for a single realization of the model inputs - was 
developed at Stanford's Center for Turbulence Research as part of the Department of 
Energy's Predictive Science Academic Alliance Program; the numerical method used 
is described in [28] and is based on an implicit, second order spatial discretization. For 
this example, we slightly modified the codes to extract of the non-zero elements of the 
matrix and right hand side used in the computation of the temperature distribution. 
With access to the linear system, we were able to apply the Galerkin method to the 
stochastic version of the problem to approximate the statistics of solution. 

4.1. Problem Set-up. The governing equation is the integral version of a two- 
dimensional steady advection-diffusion equation. We seek a scalar held <fi — <fi(x, y) 
representing the temperature defined on the domain SI that satisfies 

/ P cj> (v(s) ■ ds) = f (T(s)+T t )fo<t>-ds). (4.1) 
Jon v ' Jon v ' 

The density p is assumed constant. The velocity v is precomputed by solving the 
incompressible Reynolds averaged Navier-Stokes equations and randomly perturbed 
by three spatially varying oscillatory functions with different frequencies such that 
the divergence free constraint is satisfied; the magnitudes of the perturbations pa- 
rameterized by si, S2, and S3, respectively. We interpret the magnitudes as uniform 
random perturbations of the velocity field, which is simply an input to this model. The 
diffusion coefficient T = T(s4, S5, Sq) is similarly altered by a strictly positive, param- 
eterized, spatially varying function that models random perturbation. The turbulent 
diffusion coefficient T t is evaluated according to the Spalart-Allmarass model [28 . 
The domain J! is a channel with a series of cylinders; the computational mesh on 
the domain ft contains roughly 10,000 nodes and is shown in Figure 14.11 Inflow and 
outflow boundary conditions are prescribed in the streamwise direction, and periodic 
boundary conditions are set along the y coordinate. Specified heat flux boundary 
conditions are applied on the boundaries of the cylinders to model a cooling system. 

The goal is to compute the expectation and variance of the scalar field <f> — 
4>{x, y, s) over the domain Q with respect to the variability introduced by the param- 
eters. We use the Galerkin method to construct a polynomial approximation of <f> 
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Si 52 S3 54 S5 ■% 

3 118 5 5 

Table 4.1 

The order of univariate polynomial for each variable used in the basis set. 




Fig. 4.2. The expectation (above) and variance (below) of the temperature field <j> over the 
domain. Red corresponds to larger values and blue corresponds to smaller values. 

along the coordinates induced by the parameters s using product Legendre polynomi- 
als. To show that this method applies for an arbitrary basis set, we choose the basis 
polynomials to reflect the solution's anisotropic dependence on the parameters [6]; it 
is neither the standard full polynomial or tensor product polynomial basis. The order 
of univariate polynomial associated with each parameter is given in Table 14.11 The 
basis contains 142 multivariate polynomials. To solve the Galerkin system, we use 
Matlab's BiCGstab 34j method (since the matrix is not symmetric) preconditioned 
by the diagonal elements of the parameterized matrix evaluated at the midpoint of 
the hypercube parameter space. 

We plot the expectation and variance of <j> over the domain f2 in Figure 14.21 
which are computed in the standard way as functions of the Galerkin coefficients. 
As expected, the variance in <f> occurs in the downstream portion of the domain as a 
result of the variability in the diffusion coefficient. 

5. Summary. We have examined the system of equations arising from a spectral 
Galerkin approximation of the vector valued solution x(s) to the parameterized matrix 
equation A(s)x(s) = b(s). Such problems appear when PDE models with parameter- 
ized (or random) inputs are discretized in space, and a Galerkin projection with an 
orthogonal polynomial basis is used for approximation in the parameter space. We 
showed how the system of equations used to compute the coefficients of the Galerkin 
approximation admits a factorization once the integration is formally replaced by a 
numerical quadrature rule - a common step in practice. The factorization involves (i) 
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a matrix with orthogonal rows related to the chosen polynomial basis and numerical 
quadrature rule and (ii) a block diagonal matrix with nonzero blocks equal to A(s) 
evaluated at the quadrature points. Then matrix-vector products with the Galerkin 
system can be computed with only the action of A(s) on a vector at a point in the pa- 
rameter space; this yields a reusable interface for implementing the Galerkin method. 
The factorization also reveals bounds on the spectrum of the Galerkin matrix, and the 
Kronecker structure of the factorization gives clues to successful preconditioners. We 
tested some ideas from these prcconditioner clues on the standard test problem of an 
elliptic PDE with parameterized coefficients, and we saw that using A(s) evaluated 
at the midpoint of the parameter space was as effective as the popular mean-based 
preconditioner; the midpoint preconditioner is much easier to compute, in general. As 
a proof of concept, we applied the method to an engineering flow problem by slightly 
modifying an existing CFD code to retreive the matrix elements. 

6. Acknowledgments. This material is based upon work supported by the De- 
partment of Energy [National Nuclear Security Administration] under Award Number 
NA28614. 
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